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Abstract 


The use of genetic algorithms for minimization of differentiable 
functions that are subject to differentiable constraints is consid- 
ered. A technique is demonstrated for converting the solution of 
the necessary conditions for a constrained minimum into an un- 
constrained function minimization. This technique is extended as 
a global constrained optimization algorithm. The theory is applied 
to calculating minimum-fuel ascent control settings for an energy 
state model of an aerospace plane . 

Introduction 

Genetic algorithms for optimization (refs. 1 to 4) are nonderivative, nondescent, random- 
search procedures for functional minimization, and their algorithmic structure is based on 
biological concepts. Familiar descent-type minimization algorithms construct a sequence of 
iterations, each of which modifies the independent variable vector from the previous iteration. 
Genetic algorithms, in contrast, construct a random sequence of generations in which a 
population of codings of bounded independent variable vectors is modified according to analogs 
of biological cross breeding and mutation. Rather than finalizing the value of the independent 
variable in an iteration by satisfying a descent condition, the genetic algorithm employs a 
“survival-of-the- fittest” heuristic that assigns a greater likelihood of appearing in the subsequent 
generation to the population elements that have lower objective function values than those that 
have higher objective function values. 

A growing body of experimental evidence exists (refs. 5 to 8), supplemented by formal results 
(ref. 9), which indicates that genetic algorithms (GA’s) are reliable methods for approximately 
determining the global minimum of a function. These algorithms lack a strict descent require- 
ment, and their search operates on a population of iterates rather than on a single sequence of 
iterates. These features help prevent GA’s from becoming “stuck” at local minima. On the other 
hand, GA s do not exploit derivative information in the search. This property, coupled with the 
fact that the algorithms operate on fairly coarse codings of the independent variable vector 
(rather than on floating point numbers), tends to limit the applicability of GA’s to “rough-cut” 
analyses rather than highly accurate ones. When highly accurate solutions are required, GA’s 
can be useful to generate initial guesses for gradient or Newton algorithms. 

There have been a number of efforts in recent years to solve constrained optimization problems 
using GA’s. The most straightforward approach is to convert the constrained problem into 
an unconstrained one by adding a penalty function on the constraint violation to the cost 
function. Difficulties exist, however, which are associated with both “light” and “heavy” 
penalty weightings, just as in the case of gradient-based optimization methods. When light ^ 
penalties are employed, they generally fail to accurately enforce the constraint. When extremely 
heavy penalties are employed, that portion of the population which violates the constraints will 
have a vanishingly small probability of reproducing itself in subsequent generations. This “die- 
off” of illegal population elements results in an effectively smaller population (i.e., subsequent 
generations will have many replicates of the legal subset of the population and vanishingly few 
from the illegal subset). The resulting reduction in “genetic diversity” can adversely affect the 
performance of the algorithm. 

This paper demonstrates a GA-based approach for solving nonlinearly constrained optimiza- 
tion problems. The method, which is simple to implement and generic to structure, is applica- 
ble to problems in which the cost function and the constraints are continuously differentiable. 
Moreover, this method can be adapted to calculate the global optimizer under nonrestrictive 


assumptions. The algorithmic performance of the approach is explored in two numerical experi- 
ments. The first experiment compares the performance of the approach with a penalty function 
formulation for a simple test problem, and the second experiment extends the comparison to an 
aerospace performance optimization problem. 
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heavy penalty weight in numerical examples 
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dynamic pressure, lbf/ft 2 
space of j-dimensional real vectors 
equatorial Earth radius, ft 
reference area, ft 2 
thrust, lbf 

operator concatenating elements as vector 
vector of free parameters 
angle of attack, deg 

system of equations for stationarity of Lagrangian function 
deflection of j th control effector, deg 

parameter defining transition from light to heavy penalty weight 
fuel equivalence ratio 

number of new local minima identified in global minimization iteration 

Lagrange multiplier 

number of inequality constraints 

function returning Lagrange multiplier vector at optimum 
step-size scale factor in Newton-Raphson algorithm 
atmospheric density, slug/ft 3 
user-specified search volume 
constrained minimization function 


Subscripts: 

CG center of gravity 

( )e pertaining to elevon deflection 

( )pen pertaining to values returned by penalty function mimization 

( )t pertaining to thrust 

( )x partial derivative with respect to x 

Acronyms and Abbreviations: 

GA generic algorithm 

KT Kuhn-Tucker conditions 

max maximum 

min minimum 

NR Newton-Raphson algorithm 

Symbols with superscript stars ( )*, tildes ('), plus signs ( )+, and zeros ( )° indicate 
optimal value, active constraints, pseudo-inverse functions, and solutions returned from the 
global optimization algorithm interates, respectively. A bar above a symbol (“) indicates an 
admissible search value. 
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Problem Representation for Genetic Solution 

We treat the problem of minimizing the C 1 function c(x) subject to C 1 constraints; that is, 


x* G 7Z n is sought such that 

c(x*) < 

c(x) 

(1) 

subject to 

'x 
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/j(x*)> 0 

U € 2) 

( 3 ) 

where £ and T are the sets of indices of equality and inequality constraints, 
Identifying the active inequality constraint index set as 

respectively. 


l={j :j el,fj(x*)= 0} 

( 4 ) 

and defining 

f(x*) = vec j/ fc (x 

■*) : k efiui} 

( 5 ) 

assume that f 6 IV 1 , n < n y and 

rank {i 

} =A1 
x* J 

(6) 

If the above assumptions and equation (1) are true, then the Lagrange multipliers A exist 
(ref. 9) such that equations (2) and (3) are satisfied. 


ac(x, a) 

= 0 

x*,A* 

( 7 ) 
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A* > 0 
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( 9 ) 
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i 

II 
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Equations (2) and (3) and (7) to (9), subject to equation (6), make up a typical statement 
of the first-order necessary conditions for a constrained local minimum, or Kuhn-Tucker (K 1 j 
conditions. 

IF A* were known, a GA could be used to satisfy the KT conditions by solving for the global 
minimum (zero) of 


* (x, A*) = 53 |jCx ‘ (x ’ A * )! + ^ + 53 |min{ °’ (11) 

i=l je£ kel 

Note that, if one cast the first sum in equation (11) in the role of a cost function, then V (x, A ) 
has the structure of a typical penalty function formulation, but it has the penalty weights in the 
second two sums set to unity. The difference between the *(x,A*) and a penalty formulation 
of the constrained minimization problem is twofold. The first difference is that in the typical 
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case in which the solution is not known a priori, the optimal cost is unknown and may be 
nonzero. In this situation, it is well known that the minimum of the sum of cost and weighted 
penalty terms will vary with the selection of the penalty weighting parameters. In equation (11), 
however, all terms go to zero at x*; this situation results in the solution being invariant with 
respect to the nonunity scaling of the second two sums in equation (11). This property is 
advantageous because it eliminates ambiguity concerning the influence of penalty weightings 
on the solution. The second difference is that the “cost” in equation (11) is a measure of the 
constrained “stationarity” of the solution, rather than a direct measure of the performance. 


Because A* is not generally known a priori, consider estimating A* in equation (11) during 
execution of the GA. Define 


?( x ) = {J 'j€T, fj(x) < 0} 

as an index set of constraints which are active or violated at a given x, and 

(12) 

f(x) 

1 = vec {/ fc (x) : 

k e £ u2’(x)| 

(13) 

Now, estimate A* by fyx), where 

r i>i(x) 

(i e £) 


* 'i(x) = < 

1 i>i(x) 

(* e r(x)) 

(14) 


[o 
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and 

^( x ) = [fj( x )] + c x ( x ) 

where ( )+ denotes the pseudoinverse operator. Note that from equations (6) to (10), 

(15) 


v(x*) = 

A* 

(16) 


The use of absolute values in equation (14) is an algorithmic measure to reject constrained 
stationary points that fail to satisfy equation (8). The KT conditions are satisfied by solving 


^[x,i/(x)] = 0 (17) 

The nonsmooth equation (17) is solved in the rest of the paper by using a GA to solve the 
nonsmooth unconstrained minimization problem 

x* = argmin \P[x, i/(x)] (18) 

where x is the user-specified bounded volume over which the genetic search takes place: 

X = {x : (xj) m i n < Xj < (Xj) max (i — l,...,n)} (19) 

In this context, the GA provides a robust means to identify candidate local minima, in the sense 
of finding points that satisfy the KT conditions. Furthermore, because the minimum value of T 
is known a priori, the GA can be stopped when |\&| is “sufficiently small.” 

This latter characteristic can be useful in practice. Identification of the appropriate generation 
at which to terminate a GA minimization is an open research topic for general applications. For 
this reason, the length of the GA runs in practical studies is often set by the patience of the 
analyst or the availability of the computer resources. An objective threshold for termination can 
significantly shorten run times without loss of confidence in the solution. 
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Global Minimization Algorithm 

The use of GA’s to obtain solutions to the KT conditions, as expressed by equation (17), 
supplies the basis for a global optimization algorithm that is subject to the following two 
assumptions. Assumption 1 (i.e., smoothness) states that 


c(x) 

/i(x) 


€ C 1 xex) 


Assumption 2 states that solutions (x^,A®) of 


r(x°,A°) = 


£ x (x°,A°) 

f(x°) 


= 0 


( 20 ) 


are regular at all x° € X- These assumptions lead to the following assertion, which is proved in 
appendix A: If assumptions 1 and 2 arc true, then there are a finite number of points x £ x at 
which equation (20) is satisfied. 

The global constrained minimization procedure essentially consists of identifying all the local 
constrained inflection points {x 0 } of £ and then accepting that point which returns the lowest 
value of c(x°) as the globally minimizing solution. By this assertion, the survey will be completed 
after a finite number of identifications. 

The set x can be surveyed for the global minimizer by a penalty function-based extension of 
the approach for solving local necessary conditions. At the conclusion of a successful GA execu- 
tion of equation (18) for a given x, there will be /c > 1 roots returned, {x } {x^ , Z 1, k}, 

thus corresponding to local solutions of the necessary conditions. The problem of solving equa- 
tion (18) can then be reposed to ensure that the roots {x 0 } are excluded from the solution by 
replacing ^fx, i/(x)] with 

^'(x, {x 0 }) = 'Ffx, i'(x)] + Kcrfx, {x 0 }) (21) 

where K > 0 is a user-chosen penalty weighting term and a is a function that becomes positive 
when x is “close” to any point in the set {x 0 } but is otherwise zero. For example, 

£r= fl (xe {UF=1 £(*?)}) (22) 

( 0 (Otherwise) 

where each ^(xj*) is a user-defined volume surrounding the corresponding xj?. This approach can 
be generalized for a global survey of x by using the GA to minimize the sequence of functions 


\Pj(x) = #[x, i/(x)j -1- Kffj-i 

(j = 1,2,...) 

(23) 

{ 0"(x, Ujfc=o{ X °}fc) 

o' o 
II A 

(24) 


where {x°} fc is the set of local solutions identified when minimizing 'Fjfc-i(x). The resulting 
algorithmic structure is as follows: 

1. Set j = 1. Set costo ~ oo. 

2. Generate a random population distributed over x- 
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3. Execute a GA to minimize 'J'j (x) from equation (23). If the GA is unsuccessful in finding 
x such that 'f'j(x) « 0, go to step 8. 

4. Collect {x°}j such that ^j( x )l xe { x 0}. ~ 0. 

5. Set costj = min < costj.j, min c(x) 

1 x€{x 0 }j 

6. If cost j < costj_j, then set x* = arg costj. 

7. Set j = j + 1. Go to step 2. 

8. If j > 1, accept x® as global minimizer; if j ^ 1, abort. 

The algorithm terminates when it is no longer able to identify values of x 6 x for which ~ 0. 
It is presumed that this occurs when all the local solutions have been identified and included in 
the penalty term of tyj ’s. 

Note that the proposed algorithm would still terminate in a finite number of iterations even if 
local nonunique roots of T exist, thus invalidating the assertion. The termination occurs because 
the penalties K that exclude a growing union of finite subvolumes from On the other hand, 
the resolution of the algorithm to separately identify the closely spaced solution points depends 
on the selection of the volumes B. Recall, however, that the motivation for selecting a GA over 
a more accurate method is to robustly obtain a “rough-cut” answer. Therefore, it is not felt 
that this latter concern prevents the algorithm from having practical utility. 

The concept of using penalty functions to exclude known local minima from future iterations 
of a global optimization algorithm has been established for descent-based methods (refs. 10 
to 12). The penalty function employed in this work is similar to the “tunneling” technique 
developed in references 10 and 11. In these references, the penalty varies smoothly toward 
a huge value as x — > x/. Penalty functions of this type can lead to numerical difficulties in 
algorithms that use gradient information for defining a search direction. These difficulties are 
avoided in the present approach because of the nondescent nature of the GA. 

Numerical Experiments 

This section describes the numerical experiments that explore the performance of the GA- 
based constrained minimization procedure developed in this paper. Because the procedure 
involves more complexity than formulations in which constraints are enforced via penalty 
functions, the GA-based solutions of equation (18) are compared with the GA-based solutions 
of a “generic” penalty function problem formulation of the form 



x pen = arg min 


-( x )+ P[ x */*( x )] 


kGSUl 


where 


when k e £ and 


p[ x - /*(x)] 


B-\f k \ 

HA I 


(l/*l > e) 
(l/fel < e) 


P[ x - /fc( x )] = 


B '\fk\ (~fk > e) 

b ■ max{0, -f k } {-f k < e) 


(26a) 

(26b) 


when k € T. The parameters e > 0 and B > b > 0 are to be chosen by the user. This 
formulation allows heavy penalties that strongly violate constraints and light penalties that 
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“nearly” comply with constraints to be applied to x. Heavy penalties could be appropriate to 
reject the “artifact” local minima in equation (25) which would not approximately correspond 
to the solutions of the KT conditions for the underlying problem. The provision for lighter 
penalties within the e-defined region of light constraint penalty is intended to provide a subset 
of x in which the variation of c(x) is not dominated by the penalty terms. 

The GA used in this study was a simple GA, which was similar to that used in reference 8, 
however, it included a modification that was suggested in reference 13. In this modification, 
the best-valued population element from each generation was guaranteed survival into the next 
generation. The real- valued independent variables were coded as 8-bit binary strings such that 


x(string) — Xmin T (imax ^min) 


baseio (string) 
2 8 - 1 


(27) 


thus yielding a resolution of roughly 0.4 percent over the range of x for each problem. In both 
experiments, the string representations for the vector-valued independent variables were formed 
by concatenating the 8-bit substrings for each scalar. Key parameters affecting the performance 
of a GA (ref. 3) are population size iV pop , crossover probability P C ross, and mutation probability 
Pmutate- All runs in this study were made with Af pop = 30 and P ;ross — 0.95, as per guidelines 
from reference 3. This simple GA with the modification from reference 13, however, appears to 
benefit from a more aggressive mutation rate than the Pmutate = 0.01 that was recommended in 
reference 3. Some adjustment of this parameter was done in the experiments below. 

The first experiment (denoted example 1 in the table titles of tables 1 to 5) compares the 
performance of the generic GA penalty function approach with that of the GA solution of 
reference 11, henceforth referred to as a KT solution, for the problem 

2 2 

c{x i , x 2 ) T x 2 


f(x i,x 2 ) = XI - x 2 - 2 = 0 


with the search volume 

X = {—5 < xi < 5 ; — 5 < x 2 < 5} 

Monte Carlo experiments of 100 runs, each consisting of 100 generations, were performed 
for the KT formulation (ref. 11), and for the penalty formulation for a number of combinations 
of { P, b , e, Pmutate}- The performance results are given in tables 1 to 5. Table 1 displays the 
KT performance for Pmutate = 0.01,0.03,0.05,0.07, and tables 2 to 5, in turn, display the 
GA penalty function results for each value of Pmutate- The performance of the GA in these 
experiments was characterized by the number of runs that successfully satisfied error thresholds 
of the form err x * < k, where err x * is the Euclidean norm of the distance between the real value 
of the best population element, and the known optimal solution point x* = (1,-1). Tables 2 to 
5 also display the median number of generations necessary for the successful runs to cross the 
thresholds. 

The comparison of table 1 with tables 2 to 5 immediately reveals that the KT formulation was 
generally more successful in finding the optimum for this problem and, when successful, tended 
to find it more quickly, except in the case of Pmutate — 0-01- Note that, both for the KT and 
penalty approaches, Pmutate = 0.05 and 0.07 returned significantly better performance than the 
lower values. The case Pmutate = 0.07 did not show any significant advantage over Pmutate - 0.05 
in the penalty function runs, and it actually resulted in a small performance degradation in the 
KT experiment. If attention is restricted to the penalty results, then a closer examination of 
tables 4 and 5 suggests that all the combinations of small b (b — 10) and relatively large e 
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(e 0.1 and 0.2) tend to outperform other parameter combinations, particularly for success in 
the err x * <0.1 criterion. 


The success of these parameter combinations, which are cases 8, 9, 17, 18, 23, and 24 
of tables 4 and 5, can be loosely interpreted in terms of the processes operating in the GA. 
Although the comparatively heavy B penalties result in a population being “killed off” outside 
the region |/| < e, the high mutation rates tend to introduce enough “new genetic information” 
to prevent the population elements x from stagnating; that is, they accumulate aw T ay from the 
actual minimizing value of x simply because other population elements fall outside |/| < e and 
were eliminated by the B penalty. Similarly, the combination of larger values of e and low b 
was advantageous because it assigned more volume in the parameter space to the population 
elements that could be expected to survive through enough generations to exchange meaningful 
amounts of information through crossover operations. 


The second experiment (denoted Example 2 in the table title of table 6) extends the com- 
parison of KT and penalty function algorithm performance to a more challenging optimization 
problem: selecting optimizing altitude and control settings for an energy-state approximation 
(ref. 14) of the minimum- fuel ascent to orbit for the “Langley Accelerator” (ref. 15) aerospace 
plane concept. In this experiment, a thrust-vectoring capability is added to the model. The 
energy-state approximate solution for this problem is calculated by performing algebraic mini- 
mizations for altitude and controls along a locus of specific energies E leading to orbital injection. 
This experiment considers the algebraic minimization at a single value of E. The search variable 
ranges are 


— 1 < a < 
20000 < h < 
-20 < 6 e < 
— 20 < <5y < 
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30000 

20 

20 

1.5 


and the cost function is —dE/dm, where m is mass. At a given value of E, the nondimensional 
cost is expressed as 


c(x) 


V{E,h)I S p{v) 


cos (6t + a) 


Cpj^Se) _ 1 _ 

C T (r}) J E 


(28) 


where V = \/2 g(E — h ) and g is the gravitational acceleration, which is assumed to be constant. 
Expressions for the coefficients Isp &nd T) \ and for all constants in this problem are given in 

appendix B. Two equality constraints appear. The first constraint is a vertical acceleration 
balance 




Ci(a, 6 e ) -|- C T (rj) sin(o: + S?) 


+ m 


(“£*->})-• » 


where S is the reference area. The second constraint is a pitch moment balance 


C T (rf)x T $e)c + x C g[Cd{<*, S e ) sin o + C £ (o,<5 e )cosa]} - sin 6 T = 0 


(30) 


where q — pF 2 / 2 is the dynamic pressure and xqq and Tj are the moment arms. There is also 
an inequality constraint on dynamic pressure, which is 



9max 


> 0 


(31) 
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where g max = 2000 lbf/ft 2 . The specific energy (E = 10 5 ft) considered here approximately 
corresponds to flight at Mach 2. 

In this second example, KT and penalty function versions of the problem were again compared 
in the Monte Carlo experiments. In this case, each experiment consisted of 100 trial GA runs, 
with each run having 600 generations. Two evaluation criteria were employed to compare the 
quality of the KT and the penalty results. The first criterion was the distribution of values of 
^!(x) from equation (23); these values were returned by the populations of final-generation x’s 
from the KT and the penalty run sets. The second criterion informally quantified the usefulness 
of the GA results as initial guesses for high-accuracy Newton-type methods. This quantification 
was done by using the final-generation x’s from the GA’s as initial guesses for a restricted- 
step Newton-Raphson (NR) algorithm that solved a system of equations equivalent to the KT 
conditions for this problem: 

$i(x*) = {[/-/x(/J) + ]cx}- +fi=0 (i = l,-,5) (32) 

where /(x) is the concatenation of equations (29) to (31), with equation (31) treated as an 
equality and / r (x) = [f T (x) 00]. At the optimal solution, 

(x*) T = [1.0211, 20 385, 0.1789, 0.7182, 0.7688] 

where the treatment of equation (31) as an equality constraint is justified because the value of 
its Lagrange multiplier is positive at x*, thus taking on the value of A^ max = 6.1866. 

The NR iteration took the form 

x fc _i = X k + £ k s k s k = -[M*/t)] -1 iKzfc) (33) 

where gx{ x k) was approximated by a first-order forward difference formula, and the line search 
parameter £*. was chosen by the logic, iterated over j = 0, 1, ..., as 

(&) 0 = min{2£fc_i, 1} 

f (g[ x k i (£fc)j s fc] — 9( x k- 1)) 

iik)j = \ Abort ((&)j<£ min) 

where £ min was chosen as 10~ 7 . The algorithm was considered to have converged if the criterion 
£5 g } (x) < 1CT 4 was satisfied in 100 iterations or fewer. 

First, GA minimization of the KT formulation was considered. A Monte Carlo set of 100 
GA r uns minimizing ’Fi(x) over x was calculated. Figure 1 displays the resulting distribution 
of 'I'i values, referred to as KT errors. The median for these runs is = 0.3821, which 
corresponds to 

(4)median = [0-8353, 24 980, -1.9608, 8.8627, 0.7745] 

An examination of the final population elements from the set of trials revealed that, rather 
than clustering around a point, the 6 e and 5j components of the x’s were distributed along a 
line. This distribution is not surprising, given the correlation in the effects of 6 e and S T on 
pitching moment. There was also a significant linear trend of the deflections with a. Figure 2 
displays relationships. 

Figure 3 gives the corresponding distribution of a and tj as functions of h. The trends in 
these variables are not as strong as those seen in the relations among a, S e , and 8j. Part of the 
reason for this is that the q max constraint (eq. (31)) only weakly affects the performance for this 
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model, even though q at the maximum altitude is little more than one-half of q max . Note that 
the trends in figure 3 are explained by the reduction of q with the increase in altitude. A larger 
a is called for because of the dimunition in lift, and a slightly larger 77 is called for to balance 
the reduction in cos(< 5 r + a) and the increase in Cjj(a, 6 e ) in equation (28). 

In 100 trials, the lowest value of achieved was no better than 0.0979; this occurence can be 
explained by the properties of for this problem and by the characteristics of the crossbreeding 
and reproduction operations in GA’s. Near its minimum, ^ is least sensitive along the locus 
whose projection into (a x 8 e x <5 r )-space is depicted in figure 2. In the reproduction operation 
of the GA s, population elements away from this locus have significantly higher 'I ' 1 values and, 
therefore, are assigned a significantly lower probability of surviving into the next generation. The 
crossbreeding operation modifies pairs of population elements by swapping substrings from 
binary codings of both elements. When “near-optimum” values of are distributed along a 
surface, rather than a point, the swapping generally results in moving the modified elements 
away from the surface. These elements, in turn, lose reproduction probability and disappear 
from the population. 

To address this difficulty, 8 e and S T were transformed to tailor the search volume v to the 
behavior of such that 

6 e = £ e (a) + 6 e 8? = £ T {a) + S T (35) 

where the coefficients in 

4(a) = (ci) e a + (c 2 ) e | 

h (a) = ( c x ) T a + (C 2 )t J ^ 

were calculated by least-squares fits over the data from the Monte Carlo experiment. The values 
for the coefficients and for the residuals measure 


S( ) = 




100 

£[«( >(<*)-«( jt] 

i=I 


(37) 


were (in degrees) 


(q) e = 288.2072 
(c 2 )e = -4.2551 


(ci) T = -291.4726 
(c 2 )t = 8.6639 


as 


.6505 

s T 

= 7.9568 

II 


77 ], and the search volume was redefined 

r 0.1127 < 

a 

< 2.7832 

20 227 < 

h 

< 31 356 

-3.6505 < 

S e 

< 3.6505 

-7.9568 < 

6t 

< 7.9568 

, 0.7381 < 

V 

< 0.8194 


The bounds m y for a, h, and 77 were chosen as the sample mean values ± 1.5 times the sample 
standard deviation. The bounds on 8 e and 8 T were ±s e and ±sj, respectively. 

Again, a Monte Carlo set of 100 runs was performed. Figure 4 displays the distribution of 
l for this experiment. The median value of for this set of runs was (4M me dian = 0.0627. 
and the lowest value was (ff’iJbest = 0.0074, with the corresponding control vectors 

(xf)median = [1-0867, 23 064, 1.9988, 1.2944, 0.7699] 
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(xf)best = [1-2543, 25 770, 2.1832, 1.8149, 0.7722] 

expressed in the original x coordinates. Figure 5 gives the distribution of a, 8 e , and 8j , along 
with the boundaries of x, shown as straight lines. 

The GA-generated KT solutions were next used as initial guesses for the NR algorithm. 
Eighty-two of the NR runs converged from the first set of 100 GA solutions, based on optimization 
over the full x- Ninety-eight of the NR runs converged of the second set, based on optimization 
over the restricted x- F° r comparative purposes, the NR algorithm was executed for 100 initial 
guesses that were chosen from a uniform distribution over x- Twenty-three of the NR runs 
converged in this case. We infer that the use of a genetic algorithm to minimize T* provided an 
effective means of generating an initial guess for the NR algorithm for this example. 

Once the performance of the KT approach for this problem was established, the penalty 
scheme (eq. (25)) -was applied to the problem of minimizing equation (28) subject to equa- 
tions (29) to (31). Sequences of 100 Monte Carlo runs were performed with P mu tate = 0.05 for 
the various penalty parameter combinations from the first example. Each of the sequences of 
runs was begun from the same random number seed. Because cases 1 to 3 and 10 to 12 (all 
of which used the penalty weighting b = 10 4 ) performed poorly in the first example, they were 
eliminated from consideration in this experiment. 

Table 6 summarizes the algorithm performance for these cases. This table displays the 
minimum, first through third quartile, which is displayed as Q\,Q2,Qs, and the maximum 
values of the KT cost function 'I'i which is evaluated at the penalty function solution points. In 
addition, table 6 displays the number of successful NR convergences, denoted by N S u cc> achieved 
using the solutions as initial guesses. As was observed in the first example, cases 8, 9, 17, 18, 23, 
and 24 tended to provide better results than those of other parameter combinations, in the sense 
that Tj errors tended to be smaller and N succ tended to be larger than with other combinations. 

Cases 23 and 24 were the most successful pair in producing good initial guesses for the NR 
runs. The slight advantage seen in these latter two cases may be attributed to their employment 
of the smallest B penalty weight in the study. Because of the small B weight, the nonconstraint- 
compliant population elements are granted a somewhat higher likelihood of reproduction and 
are, thus, more likely to enhance the “genetic diversity” of the population. Nonetheless, the most 
striking characteristic of the data in table 6 is that even the best results from the penalty runs 
compare poorly with the results of the KT experiments. The former solutions are less close 
to the optimum than those of the latter, in the sense of vanishing \I>i, and they did not provide 
particularly reliable initial guesses for subsequent NR solutions. 

Figure 6 displays the detailed T i cost distribution for case 23. Figure 7 gives the distribution 
of a, 8 e , and 8j from these runs, along with the optimal solution and the case 23 solution that 
returns the lowest value of the penalty-based cost function from equation (25). Note that this 
distribution is markedly different from those in figures 2 and 4. Although the best penalty-based 
solution, 

(xj’en)best = [1-1926, 233017, 0.2354, 0.8630, 0.7591] 

is fairly close to x*, the overall trends of the penalty solutions are significantly different from 
the KT solutions. Comparing figures 2 and 4 with figure 7, note that there is a strongly linear 
trend between 8f and 8 e in all sets of solutions, but the slope of the penalty solution trends is 
opposite in sign to the KT solution trends and much different in slope magnitude. Also, the 
variation of a in the penalty function solutions is much smaller than that in the KT solutions. 
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The solution x* was tentatively verified as a global optimizer for this example by performing 
100 of the 600-generation GA runs that minimize ^2 from equation (23), using the x coordinates 
from equations (35) and (36), and 




r a - a* 
h = h* 

i ~ s e < s e 

— sj' < Sj' < 8^ -f- sj 
,77 = 77* 


This form of ^(x) was intended to deny the locus of the (£ e , 8j) pairs from the search to the 
GA. The best value of from this set of runs was 1.2715. Because this number is considerably 
higher than the worst value of ^ 1 , we infer that no other local minima were identified and, 
hence, x* is the global minimizer. 


Summary of Results 

This paper has examined the use of a simple genetic algorithm to solve minimization problems 
for differentiable functions that are subject to differentiable equality and inequality constraints. 
The first-order necessary conditions for a constrained minimum have been adapted to convert a 
given constrained minimization problem into an unconstrained minimization of a nonsmooth 
function whose minimum value is zero and whose minimization is equivalent to satisfying 
the first-order necessary conditions for the original problem. The unconstrained nonsmooth 
minimization is carried out using the genetic algorithm. 

This solution approach was exercised and compared with a penalty function formulation 
for two constrained minimization problems. In the first problem, the approach significantly 
outperformed the penalty function technique over a range of penalty function tuning parameters. 
In the second problem, the approach provided significantly more accurate solutions than the 
penalty function technique, despite numerically challenging features, such as correlated control 
variables. 


NASA Langley Research Center 
Hampton, VA 23681-0001 
March 14, 1994 
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Appendix A 

Proof of Assertion 

Assume that an infinite number of solutions x k exist in x with corresponding values of 
X k = u(x k ) from equation (15). This assumption and the fact that x is closed and bounded 
implies that {x^} will contain an accumulation point x E X- Construct a sequence of points 
x k — > x, and consider variation of ||r[x, i^(x)]|| along the lines 

x/c( a ) = (1 - (0 < a < 1) ( A1 ) 

By the extreme value theorem, for each k, there exist a k that satisfy 

K k = max ||T {x(a), i/[x(a)]}|| ( A2 ) 


As k — > oo, 




/ g«r(x 

\ dx 


x k -x 


Ofc||x fc -x|| \ dx | s ||x*-x||/ 

If K k ± 0 as k -» oo, assumption 1 of the assertion is violated. Because of assumption 1 


(A3) 


for some constant M > 0. This implies that K k -+ 0 as k -* oo, so that 

'ar(x)' 


rank 


<9x 


< n 


(A5) 


from equation (A3). Equation (A5), however, implies that there exists 0 T - [<5x r , 6X 1 ] 0 such 

that , . . 

Vr(x,A)|x i A= I /(x) 6 ' = o ( A6 ) 

which violates assumption 2. This contradiction proves the assertion. 
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Appendix B 

Smoothed Aerospace Plane Model 

This appendix describes the smooth analytical adaptation of the piecewise linear tabular 
Langley Accelerator” vehicle model described in reference 15. The following aerodynamic and 
propulsion coefficient expressions are intended for flight conditions between Mach 2 and 2.5: 

C L (a,S e ) = -0.0062 + 0.0242a - 0.000675 e + (0.896 x 10 _7 )a<5 e 

C D {a, 6 e ) = 0.0261 - 0.000206a + 0.000526a 2 - 0.000025a<5 c + 0.0000124<5 2 

C M (a, S e ) = 0.00102 - 0.0032a - 0.0002a 2 + 0.00056 e 

^SP(v) = 3713 + 1208// - 1740t? 2 

C t (v) = 0.0062 + 0.13167/ — 0.0 1 827 / 2 

where all angles are expressed in degrees. The variation of these quantities with Mach number 
has been ignored for simplicity. Figures 8 to 12 display the errors between the above analytical 
expressions and the linearly interpolated values from reference (22); the values are normalized 
by the latter and expressed as percentages. The vehicle-related constants appearing in the text 
are as follows: 

S = 3603 ft 2 
c = 80 ft 


X CG — 14.01 ft 
xt = 61.99 ft 


m = 4800 slugs 

Figure 13 displays the geometry of the vehicle. Finally, the atmospheric density model for this 
study is 


p{h) = p Q exp(-{3h) 


where p 0 = 0.002378 slug/ft 3 and 0 = 0.0000547 1/ft. 
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Table 1. Effect of Mutation Probability 
on KT Performance in Example 1 



Successes (median generations) 

^mutate 

err x * < 0.10 

err x * < 0.05 

err x * < 0.01 

0.01 

70(12) 

49(18) 

38(28) 

.03 

91(16) 

75(23) 

72(41) 

.05 

98(14) 

83(31) 

64(39) 

.07 

98(15) 

78(25) 

64(42) 


Table 2. Penalty Function GA Performance for P mu tate = 0.01 in Example 1 





Table 3. Penalty Function GA Performance for P mu tate = 0.03 in Example 1 


Case 

B 

b 


Successes (median generations) 

e 

err x * < 0.10 

err x * < 0.05 

err x * < 0.01 

i 

10 20 

10 4 

0.05 

15(51) 

8(52) 

8(52) 

2 

10 20 

10 4 

.10 

15(40) 

10(62) 

10(62) 

3 

10 20 

10 4 

.20 

18(52) 

13(73) 

13(73) 

4 

10 20 

10 2 

.05 

21(38) 

14(40) 

13(59) 

5 

10 20 

10 2 

.10 

27(39) 

18(46) 

17(56) 

6 

o 

to 

o 

10 2 

.20 

34(28) 

23(46) 

22(49) 

7 

10 20 

10 1 

.05 

29(34) 

23(62) 

23(62) 

8 

10 20 

10 1 

.10 

29(51) 

20(55) 

19(58) 

9 

10 20 

10 1 

.20 

35(35) 

27(57) 

27(57) 

10 

10 6 

10 4 

.05 

26(39) 

20(54) 

19(56) 

11 

10 6 

10 4 

.10 

27(31) 

22(53) 

21(56) 

12 

10 6 

10 4 

.20 

28(28) 

19(56) 

18(62) 

13 

10 6 

10 2 

.05 

27(59) 

19(66) 

19(66) 

14 

10 6 

10 2 

.10 

39(36) 

28(56) 

28(56) 

15 

10 6 

10 2 

.20 

37(49) 

27(62) 

27(62) 

16 

10 6 

10 1 

.05 

29(37) 

22(52) 

22(54) 

17 

10 6 

10 1 

.10 

34(54) 

27(63) 

27(67) 

18 

10 6 

10 1 

.20 

41(37) 

28(64) 

26(64) 

19 

10 4 

10 2 

.05 

26(51) 

18(68) 

18(68) 

20 

10 4 

10 2 

.10 

28(44) 

20(62) 

20(62) 

21 

10 4 

10 2 

.20 

33(39) 

23(60) 

23(60) 

22 

10 4 

10 1 

.05 

33(56) 

28(66) 

28(66) 

23 

10 4 

10 1 

.10 

37(38) 

26(66) 

26(68) 

24 

10 4 

10 1 

.20 

48(60) 

30(74) 

28(76) 
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Table 4. Penalty Function GA Performance for P muta te = 0.05 in Example 1 




Table 5. Penalty Function GA Performance for Pmutate = 0.07 in Example 1 


Case 

B 

b 


Successes (median generations) 

e 

err x * < 0.10 

err x * < 0.05 

err x * < 0.01 

i 

10 20 

10 4 

0.05 

35(52) 

18(53) 

15(60) 

2 

10 20 

10 4 

.10 

39(58) 

22(71) 

20(76) 

3 

10 20 

10 4 

.20 

38(43) 

20(52) 

20(52) 

4 

10 20 

10 2 

.05 

45(50) 

29(65) 

28(65) 

5 

10 20 

10 2 

.10 

40(46) 

27(69) 

27(69) 

6 

10 20 

10 2 

.20 

45(51) 

26(64) 

26(66) 

7 

10 20 

10 1 

.05 

47(49) 

33(46) 

29(50) 

8 

10 20 

10 1 

.10 

62(54) 

37(66) 

32(73) 

9 

10 20 

10 1 

.20 

57(38) 

38(43) 

38(48) 

10 

10 6 

10 4 

.05 

32(36) 

17(62) 

16(63) 

11 

10 6 

10 4 

.10 

41(62) 

25(75) 

25(75) 

12 

10 6 

10 4 

.20 

40(58) 

22(66) 

22(66) 

13 

10 6 

10 2 

.05 

56(48) 

28(46) 

27(49) 

14 

10 6 

10 2 

.10 

44(64) 

23(74) 

21(74) 

15 

10 6 

10 2 

.20 

49(50) 

34(66) 

33(74) 

16 

10 c 

10 1 

.05 

48(48) 

35(57) 

27(63) 

17 

10 6 

10 1 

.10 

55(38) 

24(36) 

22(40) 

18 

10 6 

10 1 

.20 

61(46) 

39(59) 

33(61) 

19 

10 4 

10 2 

.05 

39(57) 

18(67) 

16(69) 

20 

10 4 

10 2 

.10 

34(55) 

20(57) 

18(62) 

21 

10 4 

10 2 

.20 

48(52) 

27(62) 

25(65) 

22 

10 4 

10 1 

.05 

44(38) 

31(53) 

29(61) 

23 

10 4 

10 1 

.10 

57(40) 

38(58) 

36(58) 

24 

10 4 

10 1 

.20 

67(37) 

44(60) 

42(62) 
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Table 6. Summary of Penalty Function GA Performance for Example 2 
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Number of trials 


12 -| 



0 .2 .4 .6 .8 1.0 

KT error 

Figure 1. Distribution of KT error for first aerospace plane Monte Carlo experiment. 
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a. deg a, deg 5 e , deg 


Figure 2. Distribution of a, S e , and for first aerospace plane Monte Carlo experiment. 



h, ft h, ft 


Figure 3. Distribution of a, h , and rj for first aerospace plane Monte Carlo experiment. 
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CP □□ 



KT error 


Figure 4. Distribution of KT error for refined aerospace plane Monte Carlo experiment. 



a, deg a, deg 5 e > de 9 

Figure 5. Distribution of oc , 6e - &nd fc - y for refined aerospace plane Monte Carlo experiment 




KT error 


Figure 6. Distribution of KT error for penalty formulation aerospace plane Monte Carlo 
experiment. 
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Figure 7. Distribution of a, S e , and 6j for penalty function aerospace plane Monte Carlo 
experiment. 
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Figure 12. Cj- percent error in aerospace plane model. 



Figure 13. Aerospace plane geometry. Subscript CG indicates center of gravity. 
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